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Abstract 

The  paper  reviews  the  construction  of  cloaks  for  2D  acoustic  or  electromagnetic  waves  using  the 
Piola  transform,  and  shows  how  the  transform  leads  to  the  construction  of  a  quasi-optimal  test  norm  for 
the  Discontinuous  Petrov-Galerkin  (DPG)  method  with  optimal  test  functions  for  this  class  of  problems. 
Numerical  experiments  for  cylindrical  and  square  cloaks  illustrate  the  discussed  concepts  and  show  that 
the  DPG  method  is  effective  in  cloak  simulations. 
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1  Introduction 

We  begin  by  reviewing  shortly  the  concept  of  Piola  transforms  lying  behind  the  construction  of  parametric 
finite  elements  forming  the  exact  sequence  (see  [10],  p.  34). 

Piola  transforms.  Let  x  =  T(£)  be  a  smooth  transformation  of  a  domain  O  c  E3  onto  a  domain  D  C  E3. 
The  transformation  defines  a  corresponding  map  between  spaces  H 1  (Q)  and  H1(D), 

u(x )  =  u(£)  =  u(T^1(x))  =  (u  o  T^1)(x)  .  (1.1) 

The  chain  rule  implies  the  corresponding  transformation  for  gradients: 

du  du  d  A 
dxi  d£k  dxi ' 

'Email:  leszek@ices.utexas.edu.  Partially  supported  by  AFOSR  under  FA9550-09-1-0608. 
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If  we  want  to  preserve  the  exact  sequence  property,  we  need  to  transform  the  larger  space  H (curl,  Q)  using 
the  same  rule: 


In  turn,  for  the  curl  operator  we  have: 


dEk 


d 


dt 

1  dxk 


Eijk  dx  ■  ~  €ijk  dx  '  El  '  _  €ijk  +  El  eijk 


F  _  p 

OXi 


dEi  8^ 


dxj  dxk 


d2ji  _  dEi  djm  d 

dxkdxj  d^m  dxj  dxk 


(1.2) 


(1.3) 


=o 


where  e^k  denotes  the  Rizzi  symbol,  i.e., 


But, 


0  if  any  two  of  indices  ijk  arc  equal, 
e^k  =  ^  1  if  ijk  is  an  even  permutation  of  123, 

—  1  if  ijk  is  an  odd  permutation  of  123. 


dim.  dil  _  !  dxi 

€ ijk  rx  rv  J  € nml  o/.  5 

dx j  OX}*  n 


where  J  1  is  the  inverse  jacobian.  Consequently,  we  obtain  the  following  Piola  transform  for  the  curl 
operator: 


dEk  _  ,  dx 

€ ijk  ~  —  J 


dxi 


di 


dE, 


^nml 


dir , 


This  leads  to  the  transform  for  the  H (div,  fi) -conforming  fields: 


Hi  =  J 


_i  dxi 

din 


H„ 


Finally, 


dHi  _  d  /  i dx{\  ~  .dxidHkdii  _  idHk 

dXi  ~  dxi  V  dikj  k+  dik  d&  dxi  ~  '  dik 


(1-4) 


(1.5) 


(1.6) 


=o 


which  establishes  the  transformation  rule  for  the  L2-conforming  elements: 


(1.7) 


The  Piola  transforms  set  gradient  into  gradient,  curl  into  curl  and  div  into  div.  Consequently,  any  dif¬ 
ferential  operator  involving  the  grad,  curl,  and  div  operators  preserves  its  structure  at  a  possible  expense  of 
introducing  an  anisotropy  in  the  “material  data”.  In  particular,  the  Piola  transforms  preserve  the  Maxwell 
equations  and  linear  acoustics  equations. 
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2D  Maxwell  and  acoustics  equations.  Consider  the  Maxwell  equations: 


V  x  E  +  iupH  =  0, 

V  x  H  —  iueE  =  0, 


(1.8) 


where  E  and  H  are  the  electric  and  magnetic  fields,  and  e  and  //  arc  the  permittivity  and  permeability, 
respectively.  In  the  2D  transverse  electric  case,  E  =  (E\.  E^.E)  and  H  =  (0.  0.  /-/).  Using  a  clockwise 
rotation  by  90  degrees,  we  transform  the  2D  Maxwell  equations  into  the  2D  linear  acoustics  equations: 


iupp  +  div-u.  =  0 
icoeu  +  Vp  =  0, 


(1-9) 


where  pressure  p  =  H,  and  velocity  u  =  (■ u\,U2 )  =  ( E2,  —E\).  Thus  the  2D  Maxwell  and  acoustics 
problems  arc  equivalent,  and  any  result  established  for  one  case  immediately  applies  to  the  other  case  as 
well. 


Construction  of  an  invisible  cloak  for  electromagnetic  waves  is  based  on  introducing  a  singular  bijective 
map  *  =  a;(£)  that  maps  a  domain  Q  —  {P\,  where  point  E  €  U,  onto  a  domain  I)  with  a  hole,  in  such  a 
way  that  the  corresponding  inverse  map  transforms  the  boundary  of  the  hole  into  the  single  point  P.  Specific 
constructions  for  a  cylindrical  and  square  cloaks  arc  discussed  in  the  next  section. 


The  linear  acoustics  equations  with  constant  material  data  //.  e  in  domain  Q  —  P,  transform  then  into 
acoustics  equations  in  domain  D  with  anisotropic  material  data: 


iufip  +g|  =0 


(1.10) 


where 

P_  T  -  dXidXi  ! 

fl  —  /iJ,  €ln  —  6  J  .  (1.11) 

oil  din 

The  concept  is  illustrated  in  Fig.  1  showing  the  real  paid  of  pressure  corresponding  to  a  plane  wave  propa¬ 
gating  horizontally,  and  the  corresponding  composition  of  the  same  function  with  the  cylindrical  cloaking 
map  defined  in  Section  2.  The  map  maps  a  hollow  cylinder  into  the  circular  domain  in  such  a  way  that  the 
inner  circle  is  mapped  onto  the  origin.  The  rest  of  the  domain  remains  intact.  As  a  consequence  of  the 
construction,  the  plane  wave  remains  unchanged  outside  of  the  cylinder  but  it  is  modified  within  its  interior. 
The  function  on  the  right  represents  the  solution  to  2D  acoustics  or  Maxwell  equations  with  anisotropic 
material  properties  within  the  hollow  cylinder  (the  cloak)  implied  by  the  singular  map.  The  wave  does  not 
penetrate  into  the  cloaked  region  (the  inner  circle)  and  remains  unchanged  outside. 


2  Design  of  cloak  materials 

General  speaking,  there  arc  several  major  approaches  to  render  objects  invisible.  For  example,  Alu  and  En- 
gheta  [2]  proposed  to  use  plasmonic  coatings  to  cancel  the  dipolar  scattering.  But  this  technique  is  limited  to 
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I 

II 

Figure  1:  Cylindrical  cloak  in  2D.  Left:  R e(p(x)),  right:  Re(p(£))  with  p(x)  =  eluJXl  and  cloaking  map 
x  =  T(£)  from  Section  2. 

the  sub-wavelength  scale  of  the  object,  and  the  coating  depends  on  the  geometry  and  material  parameters  of 
the  object.  Milton  and  Nicorovici  [17]  discovered  that  using  a  metamaterial  coating  would  cloak  polarizable 
line  dipoles.  But  the  coating  is  affected  by  the  objects  placed  inside.  In  2006,  Leonhardt  [13]  and  Pendry 
et  al  [18]  independently  discovered  a  coordinate  transformation  mechanism  for  electromagnetic  cloaking. 
Their  mechanism  was  quite  similar  to  that  of  Greenleaf  et  «/  [11,  ?]  introduced  for  conductivity.  The  main 
idea  is  to  design  a  special  metamaterial  to  guide  electromagnetic  wave  around  the  cloaked  region.  In  May 
2006,  the  first  full  wave  numerical  simulations  on  cylindrical  cloaking  was  carried  out  by  Cummer  et  al 
[5].  A  few  months  later,  the  first  experiment  of  such  a  cloak  at  microwave  frequencies  was  successfully 
demonstrated  by  Schurig  et  al  [20]. 

Since  2006,  there  have  been  many  works  devoted  to  a  study  of  using  metamaterials  [14]  to  construct 
invisibility  cloaks  of  different  shapes  (e.g.,  [12,  16,  19,  21,  22]).  The  fundamental  idea  is  based  on  the 
principle  discussed  in  the  Introduction  in  context  of  2D  equations  -  the  Maxwell  equations  are  form  invariant 
under  coordinate  transformations. 

For  simplicity,  we  shall  consider  only  two  cloak  structures:  a  cylindrical  cloak  and  a  square  cloak. 

Cylindrical  cloak.  Following  [18],  a  cylindrical  region  r  <  Ii\  can  be  cloaked  by  a  concentric  cylindrical 
annulus  R\  <r  <  R2  through  the  following  coordinate  transformation: 

f  r'(r,6)  =  ^=^r  +  R1,  0  <r<R2l 

\  9'(r,  0)  =  9,  2  0  <  9  <  27t, 

where  the  polar  coordinate  (r,  9)  is  related  to  the  Cartesian  coordinate  (x,  y )  by 

r  =  \/x2  +  y2,  9  =  tan-1  — . 

x 

Through  some  algebraic  calculation  (details  can  be  seen  in  [15]),  we  can  obtain  the  transformation 


(2.12) 


(2.13) 
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matrix 


A  = 


dx[ 

dxj 


%^  +  ^-sin 2  6  —  ^  sin  9  cos  9  0 

U2  r 

symmetric 


^^  +  ^cos2e  o 


o 


o 


i 


(2.14) 


whose  determinant  equals 


det(d)  =  R2~Ri{R2-Ri  Ri)  =  {R2~Ri)2  r' 

K  1  R-2  y  R2  r  ’  R-2  r'  -  Ri 

Furthermore,  the  relative  permittivity  in  the  transformed  space  can  be  obtained  as 


(2.15) 


e  = 


t  t  0 

*~xx  xy  w 

J  J  n 

c yx  c yy  u 

0  o'  d 


where 


(fi2_fli)2  +  +  2l)sm20|  /det(A), 

K2  r  K-2  r 


/  _  /  _ 

^ xy  ^yx 


—  —  (2  '^>1  +  — )  sin  9  cos  9 

r  R-2  r 


/det(A), 


eyy 


R2-Ri2  .  Ri/oRi  —  Ri  .  Rl\  2a\  /a  *1 

<— R — >  +T(2-R^  +  T,COS  4/detM). 


and  d  =  1  /dct(,4).  The  permeability  y!  has  the  same  as  permittivity  d 


Square  cloak.  The  same  idea  as  above  can  be  used  to  design  a  square-shaped  cloak  with  inner  square 
width  2,S'i  and  outer  square  width  2S2.  It  can  be  seen  that  the  coordinate  transformation  [19]: 

f  x\x,y)  =  +  Si 

{  y'(x,y)  =  y(^  +  ^) 


mapped  the  right  triangle  in  the  original  space  into  the  right-subdomain  in  the  transformed  space  (see  Fig.2). 
It  is  easy  to  prove  that  the  transformation  matrix  in  this  case  is 


Ar  — 


/  S2-S1 

1  s2 

_ySi 

X 2 

0 


which  has  determinant 


S2-S1S2-S1 
det(Ar)  =  — - - ( 


S2  v  S2 

The  relative  permittivity  in  the  transformed  space  can  be  obtained  as  [15]: 


/  (dd2 


ySi  S2—S1 
$2 


0 


symmetric  (^§r)2  +  {^sf1  +  ^r)2  0  /det(Ar). 


V 


(2.16) 


(2.17) 


(2.18) 
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Figure  2:  (Left)  The  original  space;  (Right)  The  transformed  space. 


by 


Corresponding  formulas  for  the  upper,  left  and  bottom  triangles 


applying  rotation  matrix  R(9)  = 


cos  v 
sin  0 
0 


—  sin  0 
cos  6 
0 


to  the 


of  the  cloak  can  be  similarly  obtained 
right  sub-domain  with  rotation  angles 


6  =  7r/2, 7 r  and  3ir/2,  respectively. 


3  DPG  Method 

The  Discontinuous  Petrov -Gal erkin  Method  with  Optimal  Test  Functions  [7,  23]  builds  on  three  crucial 
concepts: 

•  The  idea  of  Petrov-Galerkin  method  with  optimal  test  functions  computed  on  the  fly. 

•  The  element-wise  ultraweak  variational  formulation  using  discontinuous  test  functions  that  makes  the 
determination  of  optimal  test  functions  feasible. 

•  The  choice  of  an  optimal  test  norm  that  enables  correct  “mapping  properties”  and  results  in  uniform 
stability  for  wave  propagation,  i.e.,  independent  of  the  wavenumber. 

For  reader’s  convenience,  here  we  review  quickly  these  main  points  in  context  of  linear  acoustics,  see  [23,  8] 
for  details. 


Least  squares  and  optimal  testing.  Any  well-posed  variational  problem. 


u  £  U 

b(u,  v )  =  l(v),  v  £  V 


(3.19) 
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with  two  Hilbert  spaces  U,  V,  sesquilincar  form  b(u,  v)  and  antilinear  form  l(v),  is  equivalent  to  a  linear 
operator  problem, 

Bu  =  l,  B  :  U  — »  V' ,  <  Bu,v  >v'xv=  b(u,v).  (3.20) 

What  makes  variational  problems  different  from  other  linear  problems  is  that  the  operator  B  takes  values  in 
the  dual  space  V'  equipped  with  the  dual  norm, 

\\l\\v'  =  supjplpi.  (3.21) 

v^O  IlGlV 

Given  a  finite-dimensional  subspace  Uh  C  U ,  the  least  squares  method  seeks  the  minimizer  of  the  quadratic 
minimum  residual  problem, 

\\\Buh  ~  l\\v'  min  •  (3.22) 

z  uhe  uh 

Recalling  the  notion  of  the  Riesz  operator  for  test  space  V, 

Ry  '■  V  — >  V' ,  <  Ryv,5v  >=  (v,Sv)y  (3.23) 

and  the  fact  that  Ry  is  an  isometry,  we  can  reformulate  the  least  squares  problem  as, 

t§Rv1(Buh  ~  l)\\v  ->•  mm  •  (3.24) 

z  uheuh 

The  corresponding  necessary  and  equivalent3  condition  for  the  minimizer  Uh  takes  form  of  the  linear  varia¬ 
tional  equation: 

{Ry\Buh  -  l),  RylB5uh)y  =  0,  Suh  <E  Uh  (3.25) 

or,  recalling  the  definition  of  the  Riesz  operator, 

<  (Buh- l),Ry1B5uh>=  0,  5uheUh.  (3.26) 

Defining  trial-to-test  operator  with  the  corresponding  image  identified  as  the  optimal  test  space  14, 

Uh  3  Suh  — *■  vh  :=  R^BSuh.  G  Vh  C  V.  (3.27) 

we  can  reinterpret  condition  (3.26)  as  a  Petrov-Galerkin  method, 

<  (Buh- l),vh>=  0,  Rvvh  =  BSuh ,  (3.28) 


or,  returning  to  variational  notation. 


r  Uh  ^  uh 

{  b(uh,vh)  =  l(vh),  Vh  €  Vh, 

’Under  the  assumption  of  the  well-posedness,  the  least  squares  quadratic  functional  is  positive  definite. 


(3.29) 
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where  the  test  functions  arc  obtained  by  solving  an  auxiliary  variational  problem  (inverting  the  Riesz  oper¬ 
ator  Ry), 


(  vheVh 

\  (vh,Sv)v  =  b(6uh,  Sv),  SveVh. 

The  main  points  of  the  Petrov-Galerkin  method  arc: 


(3.30) 


•  The  stiffness  matrix  is  always  hermitian  and  positive-definite. 

•  The  method  delivers  the  best  approximation  error  (BAE)  in  the  “energy  norm”: 

1 1 it 1 1 £  :=  1 1 .Bit 1 1 =  sup  Ly-  ’  — .  (3.31) 

vev  IMIv 

•  The  energy  norm  of  the  discretization  error  u  —  Uh  equals  the  residual  and  can  be  computed  without 
knowing  the  exact  solution, 

\\u  —  Uh\\E  =  || Bu  —  Buh\\v  =  IK  ~  Buh\\v  =  H-Ry1^  —  Buh)\\y.  (3.32) 

Hence,  there  is  no  need  for  a  separate  a-posteriori  error  estimation,  the  method  comes  with  a  built-in 
error  estimator4. 


In  fact,  we  do  not  have  a  single  method  but  a  family  of  such,  any  choice  of  test  norm  ||u||y  results  in  a 
different  Petrov-Galerkin  approximation.  An  obvious  question  to  ask  is  how  to  select  an  optimal  test  norm. 
A  possible  answer  comes  from  the  Banach  Closed  Range  Theorem.  Under  the  assumption  that  the  dual 
operator  is  injective,  if  we  define  the  test  norm  as: 


\v\\v  ■=  sup 
ueu 


\b{u,v)\ 


Mu 


(3.33) 


the  corresponding  energy  norm  coincides  with  the  original  norm  in  U, 


Me  =  INI  u- 


(3.34) 


We  deliver  the  best  approximation  error  in  the  norm  of  our  choice. 


Ultraweak  variational  formulation  and  the  DPG  variational  formulation.  We  review  now  the  main 
ideas  of  the  DPG  method  of  C.  Bottaso,  P.  Pausin,  S.  Micheletti,  and  R.  Sacco  [3,  4]. 


We  solve  the  linear  acoustics  equations: 


iuu  +  Vp  =  0 
iup  +  divw  =  0 


(3.35) 


4Note  the  connection  with  implicit  a-posteriori  error  estimation  techniques  aiming  at  element-wise  approximation  of  inverse 
Riesz  operator  [1], 


in  a  bounded  domian  17,  accompanied  by,  e.g.  hard  boundary  condition: 

un  =  u  •  n  =  g  on  T  =  <917.  (3.36) 

We  denote  17/,  for  a  disjoint  partition  of  17  into  open  elements  K  such  as  17  =  U KenhK,  see  Fig.  3. 


Figure  3:  A  FE  mesh  with  elements  K,  edges  e,  skeleton  F/,  and  internal  skeleton  r9. 


For  any  element  K,  multiplying  the  equations  with  test  functions  v  e  H(dlv.  K).  q  E  //’  (K)  and 
integrating  over  the  element  K,  we  obtain 


iui  u  ■  v  +  Vp  ■  v  =0 
JK  J  K 

iu>  /  pq+  divw  q  =0. 
Jk  Ji< 


We  integrate  by  parts  (relax)  both5  equations: 


iu  u  ■  v  —  p  ■  div-u  +  /  pvn  =  0 
Jk  Jk  JdK 

iu  /  pq—  /  u-S7q+  /  un  q  sgn(n)  =  0, 

,  Jk  Jk  Jai< 


where  u„  =  u  ■  nP  and 


sgn(n)  = 


1  if  n  =  ne 
—1  if  n  =  —nP. 


(3.37) 


(3.38) 


(3.39) 


Finally,  contrary  to  the  concept  of  numerical  flux,  we  declare  both  traces  and  fluxes  to  be  independent 
unknowns: 


iu>  u  ■  v  —  /  p  ■  divy  +  /  pvn  =  0 
Jk  Jk  JdK 

iuj  /  pq—  u-Vq+  /  un  q  sgn(n)  =  0. 
,  Jk  Jk  Jdi< 


(3.40) 


5Hence  the  name  of  the  ultra-weak  variational  formulation. 
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We  use  BCs  to  eliminate  known  fluxes 


ico  u  ■  v  —  p  ■  divu  +  /  pvn  =  0 
Jk  Jk  JdK 

iu  pq-  u-Vq+  unqsgn(n)  =  g  q, 

Jk  Jk  JdK- r  JdK  nr 

sum  up  over  all  elements,  and  introduce  the  standard  DG  notation: 

iuj(u,  v)nh  -  (p,  di\v)nh+  <  p,  vn  >rh  =  0 
iuj(p,  q)nh  -  (u,  Wq)nh  +  <  un,  q  >ro  =<  g,  q  >r, 


(3.41) 


(3.42) 


where  we  denote  (r,  s)nh  :=  (r-  s)  K,  which  reflects  the  element  by  element  calculations.  Further¬ 

more,  the  precise  functional  setting  is  as  follows: 


Th 

Th 

tf1/2(  rh) 

h\\m/2(rh) 

H~1/2(T0h) 


:=  U  KdK  (skeleton), 

:=  Tft  —  (Kl  (internal  skeleton), 

:=  {or|rh.  :  Q  F  H1^)}  with  the  minimum  extension  norm: 

:=  inf{||Q||tfi  :  Q\rh  =  ?}, 

:=  { vn | rh  :  v  E  H o(div,  O) }  with  the  minimum  extension  norm: 
:=  inf{||v||_jjo(diVjQ)  :  v  ■  n|ro  =  an}. 


We  have  two  group  variables: 

Solution  U  =  (u,p,un,p): 

ui,u2]p  6  L2(Qh), 
un  E 

p  e  iF1/2(rfe), 

and  test  function  V  =  (v,q): 

v  E  H(drv,Uh), 
qeH 


The  sesquilinear  form  can  be  written  as: 


(3.43) 


(3.44) 


(3.45) 


b(U,  V )  =  -(u,  iuv  +  Wq)nh  -  ip,  iuq  +  divv)nh 
+  <un,q  >ro  +  <  p,  vn  >rh  ■ 

There  are  two  main  punchlines: 


(3.46) 


•  Local  invertibility  of  Riesz  operator.  Due  to  the  use  of  “broken”  Sobolev  spaces  (discontinuous  test 
functions),  the  Riesz  operator  is  inverted  elementwise.  Given  any  trial  shape  functions,  we  compute 
the  corresponding  optimal  test  functions  on  the  fly. 
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•  Approximate  optimal  test  functions.  The  locally  determined  optimal  test  functions  still  need  to 
be  approximated.  This  is  done  using  standard  Bubnov-Galerkin  method  and  an  enriched  space.  If 
polynomials  of  order  p  are  used  to  approximate  the  unknown  velocity  and  pressure,  the  approximate 
optimal  test  functions  are  determined  using  polynomials  of  order:  p  +  A p,  where  A p  >  1 . 


Requesting  the  L2-norm  for  velocity  u  and  pressure  p,  and  some  unknown  (at  this  point  yet)  norms  to  define 
the  minimum  extension  norms  for  traces  and  fluxes, 


\\{u,p,un,p)\\lr  =  ||u|||2  +  \\p\\2L2  +  H«n||?  +  ||p||?, 


(3-47) 


we  can  compute  the  optimal  test  norm  (3.33): 


ll(«.9)ll  opt  =  +  Vq\\nh  +  \\iu>q  +  divvWlb 


+  sup. 


Un,P 


,q>  +  <p,Vn>  | 

dl?  +  l|P?)1/2  ' 


(3.48) 


The  terms  involving  traces  are  unfortunately  non-local6,  which  makes  the  optimal  test  norm  unsuitable  for 
computations. 

The  following  quasi-optimal  test  norm  was  proposed  in  [23,  8]: 


\\(v,q)\\ Ipt  =  \\iu>v  +  Vq\\lh  +  \\iuq  +  divw||&h  +  |M|2  +  ll<?f ,  (3-49) 

where  ||  •  ||  denotes  the  L2-norm  on  O.  With  the  same  norm  used  to  define  the  minimum  energy  extension 
norms  for  traces  and  fluxes,  it  was  shown7  in  [8]  that  the  quasi-optimal  test  norm  is  equivalent  to  the  optimal 
test  norm  uniformly  in  wave  number  k  and  mesh  parameters:  element  size  h  and  polynomial  order  p. 

Consequently,  we  have  the  robust  stability  in  the  desired  norm: 

(||«  -  uh\\2  +  \\p~Ph\\2  +  II Un  unh || 2  T  \\p-Ph\\2)2 

<  ||  {u,p,un,p)  -  (uh,ph,Unth,ph)\\E 

=  BAE  of  (u,p,  un,p)  in  the  energy  norm 

<  BAE  of  (u,p,  un,p)  in  the  desired  norm. 


The  ultimate,  pollution-free  DPG  method  for  wave  propagation  problems.  Despite  the  uniform  stabil¬ 
ity,  the  result  above  does  not  prove  that  DPG  is  a  pollution  free  method  in  multidimensions.  This  is  because 
the  best  approximation  error  term  includes  the  best  approximation  error  for  traces  and  fluxes  which  is  not 
pollution  free.  In  ID,  fluxes  and  traces  arc  just  numbers,  so  their  best  approximation  error  (in  any  norm)  is 
always  zero. 

6For  regular  functions,  they  can  be  interpreted  as  jumps  in  vn  and  q  on  I\. 

7Under  standard  technical  assumptions. 
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The  method  proposed  and  analyzed  in  [9]  is  a  slight  variation  of  the  presented  method,  and  it  is  obtained 
formally  by  rescaling  the  L2-term  in  (3.49)  with  a  small  constant  a, 

ll(w  9) II opt  =  ¥“V  +  Vg||nh  +  \\iuq  +  divw||&h  +  a  (||v||2  +  ||<?||2)  (3.51) 

and  passing  with  a  — >  0.  The  corresponding  method  is  pollution  free,  see  [9]  for  the  proof.  In  practice,  we 
use  the  test  norm  above  with  a  small  value  of  a  limited  by  round  off  error  only. 


4  DPG  Method  for  Cloaking  Problems 


We  consider  the  2D  acoustics/Maxwell  problem  discussed  in  the  previous  section: 


iuj  u  ■  v  —  p  ■  div-u  +  /  pvn  =  0 

J  J  J  dQh 


iw  pq-  u-Vq+  un  q 

J  Qfa  J  Qh  J  d^h — r 


sgn(n)  =  / 

Jd 


(4.52) 


gq- 


8Qh  nr 


The  quasi-optimal  test  inner  product  discussed  in  Section  3  has  the  form: 


/  (iujpq+drvv)(iujpSq  +  divdv)dx+  /  (iuev+Vq)(iueSv  +  V5q)dx+a  I  /  pqSq  dx  +  /  ev5v  dx 
j  j  \j  Qji  j  Qh 

(4.53) 

As  usual,  the  differential  operators  are  understood  elementwise  (as  indicated  by  integration  over  Flf).  Note 
that  we  have  used  e  and  p  to  define  the  L2  weighted  norms. 

The  cloaking  stretching  x  =  x(£)  results  in  the  test  inner  product  in  the  £  plane, 


dik  d£i 


dq 


ddq^ 


j  (iujpq  +  divv)(iup6q  +  div5v)  J  ldf  +  J  ^ J (iujeknvn  +  ^-)(iueim5vm  +  d£ 


where 


+a  (  pqSq  d£  +  I  enivn6vi  d£  )  , 

Ifih  J  tih 


„  dxi  dxi  x 
p  =  pJ,  eni  =  e—  —  J  , 


(4.54) 

(4.55) 


d^n  <96 

and  also  the  metric  in  the  second  term  has  changed  with  the  new  metric  tensor  defined  by: 

djk  d£i  j 

dxi  dxi 

Thus,  if  we  believe  in  the  quasi-optimal  test  norm  before  the  stretching,  the  formula  above  defines  the  right 
test  norm  to  be  used  for  the  cloaking  problem. 


(4.56) 


4.1  Basic  Numerical  Experiments 

In  this  section,  we  will  use  our  DPG  method  to  solve  two  cloaking  problems. 
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We  use  the  following  data: 


Cylindrical  cloak 

wave  number: 
domain: 

cloak  dimensions: 
material  constants: 


k  =  20.958 

n  =  (-1.5, 1.5)2 

R\  =  0.15,  i?2  =  0.3 
e  =  l,ii  =  l. 


Consistent  with  discussion  in  Section  2,  we  use  the  hard  boundary  condition  along  the  inner  radius  of  the 
cloak.  Along  the  outer  boundary  of  the  domain,  we  use  impedance  boundary  condition, 

p-un  =  g  (4.57) 

with  data  g  driving  the  problem  and  determined  using  the  very  plane  wave  that  we  simulate.  Impedance 
boundary  condition  (4.57)  defines  a  relation  between  traces  p  and  fluxes  un  on  the  exterior  boundary  YfVXt. 
We  have  implemented  it  using  Lagrange  multipliers.  More  precisely,  we  enforce  the  following  constraints 
on  the  unknown  traces  and  fluxes, 

[  ip-  un)(t>  =  f  g4>  (4-58) 

^/pext  J  Y'ext 

for  every  test  function  <i>  coming  from  the  flux  space.  With  equal  order  of  approximation  for  traces  and 
fluxes  being  used  in  the  simulations,  the  condition  is  equivalent  to  a  pointwise  condition, 

p-un  =  Phg  on  Text,  (4.59) 

where  P\,  is  the  L2-projection  on  the  flux  space  on  Yext. 

We  use  the  finite  element  mesh  displayed  in  Fig.  4  with  roughly  four  bilinear  elements  per  wavelength. 
The  color  bar  on  the  right-hand  side  represents  order  of  H 1  -conforming  elements,  in  this  case  a  uniform 
order  p  =  2.  Consequently,  the  order  of  ^-conforming  elements  8  used  to  discretize  “field  variables”: 
pressure  p  and  velocity  components  u\ .  u-i  is  equal  one  (bilinear  elements). 

The  shape  of  finite  elements  reflects  trans finite  parametrizations  used  to  model  the  geometry  (see  [6] 
for  details).  Traces  arc  discretized  with  continuous  (and  hence  /7  l  2-contorming)  elements,  and  fluxes  arc 
approximated  with  quadratic  discontinues  elements.  The  optimal  test  functions  arc  determined  using  the 
enriched  space  with  A p  =  2.  More  precisely, 

q,  5q  €  Vp+Ap  <g>  Vp+Ap ,  v,5v£  ( Vp+Ap  ®  Vp+Ap~ X)  X  (pP+AP“1  ®  Vp+Ap ) 

with p  =  2  and  A p  =  2.  We  have  used  a  =  l.d  —  5.  Using  a  =  1  gives  only  slightly  bigger  L2-errors. 

The  real  paid  of  the  pressure,  shown  in  Fig.  5,  differs  very  little  (in  an  eye-ball  norm...)  from  the  exact 
solution  presented  earlier  in  Fig.  1.  The  only  visible  distorsions  arc  attributed  to  the  use  of  non-affine 
elements  and  resulting  non-polynomial  shape  functions  that  affect  the  best  approximation  error. 

8The  code  supports  the  whole  exact  sequence,  i.e.  H1-,! 7(curl)-,  H (div)-  and  L2-conforming  elements,  with  the  order  dictated 
by  the  exact  sequence  for  Nedelec’s  triangles  and  quads  of  the  first  type. 
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Figure  4:  Cylindrical  cloak  in  2D:  The  FE  mesh. 


Square  cloak.  We 

use  the 

same  data  as  for  the  cylindrical  cloak: 

wave  number: 

k  = 

20.958 

domain: 

n  = 

(-1.5, 1.5)2 

cloak  dimensions: 

Si  -- 

=  0.15,  S2  =  0.3 

material  constants: 

e  = 

l,/r  =  l. 

We  also  use  the  same  mesh  with  roughly  four  bilinear  elements  per  wavelength.  Fig. 6  presents  exact  and 
approximate  solutions  for  the  problem.  The  relative  L2 -error  for  the  field  variables  (u.  p)  is  23.77  percent. 

4.2  Experiments  with  h- Adaptivity 

To  reduce  the  relative  errors,  we  further  developed  the  h-adaptivity  for  our  DPG  method. 

We  first  tested  the  “square  cloak”  problem  to  see  how  this  works.  Our  initial  mesh  has  rougly  four 
bilinear  elements  per  wavelength,  and  our  h-adaptivity  has  produced  essentially  uniform  refinements.  The 
results  show  0(h2)  convergence  rate  and  the  ratio  of  L2  and  energy  errors  stays  bounded,  which  are  consis¬ 
tent  with  our  recent  theoretical  result  [9].  Detailed  results  are  presented  in  Table  1,  where  DOFs  denote  the 
total  Degrees  of  Freedoms. 

Similar  results  have  been  observed  for  the  “cylindrical  cloak”  problem,  details  are  presented  in  Table  2. 
The  obtained  optimal  mesh  after  the  fourth  refinement  is  shown  in  Fig.  7.  We  also  present  the  numerical  and 
analytical  solutions,  and  the  error  in  Fig.  7. 
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Figure  5:  Cylindrical  cloak  in  2D:  The  real  part  of  the  computed  pressure,  range:  (-1.22,1.31). 


Table  1:  The  “square  cloak”:  Flistory  of  h-refinements. 


Mesh 

DOFs 

Energy  error 

L2  error 

L2  relative  error  (in  %) 

L2/energy  error 

1 

2304 

1.0922923104 

0.3963140791 

9.2788672579 

0.3628278578 

2 

11064 

0.3425376653 

0.1002458536 

2.3470467350 

0.2926564398 

3 

45064 

9.1168360988E-02 

2.54643 11229E-02 

0.5961935165 

0.2793108371 

4 

178512 

2.33 170277 18E-02 

6.6201235622E-03 

0.1549963284 

0.2839179865 

4.3  Experiments  with  Pollution  Errors 

Finally,  to  demonstrate  that  our  DPG  method  exhibits  no  dispersion  error  as  proved  in  [9] ,  we  solved  our 
cloak  problems  with  different  wave  numbers  and  different  mesh  sizes  (but  keep  kh  fixed).  Obtained  numer¬ 
ical  results  are  presented  in  Table  3,  which  clearly  illustrates  lack  of  pollution  errors. 


Table  2:  The  “cylindrical  cloak”:  History  of  h-refinements. 


Mesh 

DOFs 

Energy  error 

L2  eiTor 

L2  relative  error  (in  %) 

L2/energy  error 

1 

2304 

1.4105574668 

0.4205521701 

9.8782505902 

0.2981460734 

2 

9348 

0.4918732768 

0.1359146653 

3.1924660088 

0.2763204910 

3 

12572 

0.4434205623 

0.1076659299 

2.5289206723 

0.2428077068 

4 

14282 

0.3988570480 

0.1057855040 

2.4848211457 

0.2652215990 
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Figure  6:  Cylindrical  cloak  in  2D:  The  real  part  of  computed  pressure.  Left:  the  exact  solution,  range: 
(-1,1);  Right:  the  approximate  one,  range:  (-1.24,1.29). 


Table  3:  Experiments  with  pollution  errors  for  both  cloaking  problems. 


square  cloak 

cylindrical  cloak 

Wavenumbers 

DOFs 

L2  relative  error  (in  %) 

L2  relative  error  (in  %) 

0.52396E+01 

4800 

9.27887 

9.87825 

0. 10479E+02 

18960 

9.32685 

9.68599 

0.20958E+02 

75360 

9.39416 

10.07621 

0.419 17E+02 

300480 

9.49858 

10.19268 

5  Conclusions 


In  this  paper  we  extend  our  recently  developed  Discontinuous  Galerkin  Method  (DPG)  to  the  cloaking 
problems.  Extensive  numerical  experiments  demonstrate  that  the  DPG  method  is  indeed  free  of  pollution 
eiTor  and  quite  efficient  in  solving  this  type  of  problems. 
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Figure  7:  Cylindrical  cloak  after  4  refinements.  Top  Left:  The  mesh;  Top  Right:  The  real  part  of  computed 
pressure,  range:  (-1.01,  1.05);  Bottom  Left:  the  exact  solution,  range:  (-1,1);  Bottom  Right:  the  error 
approximate  one,  range:  (-0.057,0.037). 
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